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y— i Several ways to accelerate the solution of 2D/3D linear min-max problems in n constraints are discussed. We also 
present an algorithm for solving such problems in the 2D case, which is superior to Cgal's linear programming solver, 
both in performance and in stability. 

^ 1 Purpose 

SO This work is focused on several ways to accelerate the solution of 2D/3D linear min-max problems in n constraints. We 

^N also present an algorithm for solving such problems in the 2D case, which is superior to Cgal's linear programming solver, 

! , both in performance and in stability. 

Problem 1. Linear Min-Max problem, also known as linear optimization 
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where M — \a\b\c\. 

This problem can be re-written as a Linear Programming of the following form. 
Problem 2. 

minimize.,, v t t 



in 

C*^ s.t a\x + b\y + c\ < t 

(N 

>- a n x + b n y + c n < t. 
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It is known for several decades ( |Megiddo ( 1984| ) that general Linear Programming problems (and thus, also linear min- 



^ max problems) can be solved in O (n) time when the dimension is constant. Therefore, there is only hope to demonstrate 
a constant factor acceleration; however, in real applications, this can be valuable. 



2 The Hough Transform 



We begin with a quick overview of a natural extension to the Hough Transform Hough (1962) to 3D. We note that the 
following theorems, although proved for 3D, hold equally well in 2D; the proofs are identical if we rename the z coordinate 
into y. 

Given a point p — (a, b, c) in K 3 , we define its dual plane as % (p) — I (a;, y) — ax + by — c, and given a plane 
7r (x, y) — ax + by + c, we define its dual point as % (it) = (a, 6, — c). The usefulness of these definitions is highlighted in 
the following lemmas. 

Lemma 3. A point p = (xo, yo, Zo) is above a plane tt (x, y) = ax + by + c iff the plane T-L (p) is below the point "H (it). 
Moreover, p is on n iff % (it) is on % (p). 
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Figure 1: Hough Transform of three points and a line (plane in 3D). Lines and points, and their duals, are related by their 
color. 



Proof. We know that it (xo,yo) = x oo- + Hob + c < z o] by definition of the Hough transform, \H (p)] \ T~L (tt) x ,T~L (^) y J = 
xqH (ir) x + i/qH (it) — zq which equals xqcl + yob — zq. Because xqol + yob + c < zq, we must have xoa + yob — zo < —c, 
concluding that \H (p)] {% (it) x ,H {^y^j < H (^) z - A very similar argument can be used to show that a point is on a 
plane iff its dual is on its dual. □ 

Lemma 4. The upper envelope of a set of planes corresponds to the lower convex hull of the planes' dual points. 

Proof. By definition, a point p on the upper envelope of a set of planes {fti} is above all of them (except several on which 
it lies), and by Lemma [3] we have that the plane % (p) is below all the points % (tti) (except several which it touches). 
The converse is also true - every plane it that is part of the lower convex hull is lower than all points pi (except several 
which it touches), therefore H (it) must be above all the planes H(pi). □ 

We will henceforth denote the convex hull of a set of points P by CH (P), and its lower convex hull by CH (P). Figure 
Q] demonstrates the above lemmas in the 2D case. 



3 A Linear Min-Max Problem in the Hough space 

Solving the linear program is equivalent to finding the lowest point of the upper envelope of the set of planes defined by 

Zk (ar, y) = a k x + b k y + c k=l,2,...,n. 

We saw that the upper envelope corresponds to the lower convex hull of the set of points DP = {(a*;, b k , — Cfe)}? =1 . It 
is now obvious that a solution to Problem [2j which is a point (x : y,t) on the upper envelope, corresponds to a plane in 
the dual space, which is defined by a face of CH (DP). 

Consider an optimum for the target function of the linear problem, namely t opt , and define the following plane which 
is parallel to the x — y plane: n (x, y) = • x + • y + 1 • t opt . The dual to this plane, % (it), is a point on the z axis. Because 
the optimal solution p to Problem [2] is a point that is on n, its dual T-L (p) is a plane on which the point T-L (tt) must be. 
This means that the dual to the plane that has the highest intersection with the z — axis (and is part of CH (DP)) is 
the solution to the linear program. Moreover, it is obvious that if there is a face of the lower convex hull that intersects 
the z axis, it has the highest intersecting plane. If there is no such face, either all of the points have a positive a k , or all 
have a negative a k , which means the problem is unbounded. 

The last two statements suggest that solving our linear program is equivalent to finding the face of the lower convex 
hull that intersects the z axis. Unfortunately, we are unaware of an O (n) method to do so. 

Figure [2] demonstrates the above lemmas in the 2D case. 
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Figure 2: Left: a linear min-max problem, and its solution (red point). Right: duals of planes defining the problem, and 
the dual to the solution (red line). 



4 The 2D Case: Solving Problem [2] in the Hough plane 

Although the last section ends in a pessimistic note, this section provides an algorithm which in practice, solves 2D linear 



min-max problems much faster than Cgal's SOLVE_linear_program Fischer et al. (2011). We were unable, at the 



time of writing this work, to provide a complexity proof; however, experiments support the conjecture it is linear in the 
number of constraints. See Figure [3] 

To clarify: we remove the y— coordinate from our problem formulation and rename z into y, to obtain the following 
class of problems: 



minimize^ t 

s.t a%x + bi < t 

a n x + b n <t. 



The algorithm is as follows. 
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Figure 3: CPU time in milliseconds to solve 10,000 random problems with n constraints (our algorithm) 



Algorithm 1 Solving in Hough plane (2D) 



Input: A set of n constraints: a^rr + bi < t, i G [n] 
Output: Optimal primal point (x,t) 

1 . Transform constraints into a set of points DP 

2. Partition DP into L (points with negative x— coordinate) and R (positive) 

3. Pick some point po £ L 

4. Repeat until no change: 

(a) Given pi in L (R), find a point Pi+i in R (L) with the largest clockwise (counter-clockwise) turn from pi. 



5. Compute the line that intersects the last two points p^ and Pk+i- £(x) = mx — mp^ + p)? 1 where 



Ay) 



6. Compute the dual to this plane: (m,mp^ — p^^j 

7. Return either (6) or — oo. 



First, we should prove the algorithm takes a finite number of steps. 
Claim 5. Step (4) in Algorithm [I] terminates. 

Proof. Suppose we are at step i + 1, and the last two points are pi and Pi+i- Also assume without loss of generality 
that pi £ L and therefore Pi+\ £ R. Next, either step (4a) terminates (we're done) or a new point Pi+2 € L is selected. 
The criterion for the selection is that Pi+2 constitutes a counter-clockwise turn from pi , which is equivalent to to the line 
(Pi+2)Pi+i) having a larger inclination than the line (j>i,Pi+i). However, because both lines intersect Pi+i, and the second 
one has a larger inclination, its intersection with the y axis is smaller than that of the first one. This means that the 
y— intersection of any (pi,Pi+i) is smaller than the y— intersection of any line (pj,pj + \) for all i > j. Because the number 
of points is finite, the minimum of y— intersections of lines defined by any pair of points in DP is finite, therefore step 
(4a) terminates. □ 

We are left with showing that when step (4a) terminates, the points ph and pk+i are indeed part of CH (DP); the fact 
that the segment (j>k,Pk+i) intersects the y axis is obvious. 
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Figure 4: Two applications of step 4a 



Proof. Assume without loss of generality that pt £ L and therefore pk+i £ R- By step (4a) we know that no point in R 
has a clockwise turn from the line (pk^Pk+i), or equivalently, that all the points in R are above this line. Similarly, no 
point in L is above this same line for symmetric reasons. This means this line supports DP (from below), which means 
it must be on its lower convex hull. □ 

Figure [4] demonstrates the algorithm. 



4.1 Experimental Verification 



We pit our algorithm against Cgal's linear programming solver Fischer et al. (20111. The constraints were drawn from 
a 2D Gaussian random variable X ~ [JV (E = 0, a 2 = 10)] . As is evident in Figure 5 which depicts the ratio between 
running times of the two solvers, versus the number of constraints n, the new algorithm is about 10 times faster than 
Cgal's. 



4.1.1 Implementation notes 

Our solver only uses addition and multiplication and is highly-parallelizable. No divisions are used, which results in a 
faster and more stable solver. 

Cgal's solver requires the use of so-cold exact types, such as rational numbers or arbitrary-precision floating point 
numbers; the use of such types is extremely slow at this time, because of the way Cgal uses them. Therefore, the solver 
was tricked to using simple machine double-precision floating point numbers. This resulted in numerical failure when the 
number of constraints n was greater than w 10 5 . Our solver, however, is working even with this number of constraints. 

Moreover, the only numerically sensitive step in our algorithm is in Step (4a), which involves deciding if three points 

? 

define a clockwise turn. This is done by computing (pi — po) A (j>i — po) > 0, which can be done by setting Pi = Pi — po 
and then deciding if 

Pi ■ p\ is greater than p\ ■ p\ , 

or equivalently, 

~yC J-.X 

— =r is greater than 
Pi Pi 

This comparison can be performed in exact by converting each fraction to a continued fraction, which is still much faster 
than using an arbitrary-precision floating point number. 

A computer using Core 2 Quad (Q9400) @ 2.66 GHz and 6GB RAM was used for the benchmarking. 10,000 random 
problems were solved by each of the solvers, and their results compared. 
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Figure 5: Ratio of time to solve 10,000 problems, vs. problem size (n) 

5 The 3D Case: Discarding of Constraints 

To again recede to a pessimistic note, we were unable to extend the 2D algorithm into 3D. However, modifying the problem 
at hand allows us to quickly and safely discard of some of the constraints. The modified problem has two more constraints: 
the solution point (x, y,t) must have its coordinates x and y between and 1. 



Problem 6. 



minimize^. t t 

s.t a%x + biy + c\ < t 



a n x + b n y + c n <t 
x G [0, 1] 

y e [o, l] . 



This problem is indeed very specific, but is in the core of the GMDS algorithm Bronstein et al. (2006 1 in the norm, 
where x and y represent barycentric coordinates inside a triangle. 

To reiterate, solving Problem |6]is equivalent to finding the face of the 3D lower convex hull of the points which are dual 
to the planes defined by the constraints, which intersects the z axis. In addition, if this (only) plane z (x, y) = ax + by + z 
has either a £ [0, 1] and/or b ^ [0, 1], the solution must be on the boundary of the feasible set. 

In other words, it is possible to discard of all of the points in DP which support planes with either a ^ [0, 1] and/or 
b £ [0, 1] and not change the solution, provided that the boundaries of the feasible set are checked. We note that a 
corresponds to the inclination in the x direction, and b to the inclination in the y direction. 

To illustrate the points which can be discarded, consider a simplification of the problem to 2D. Figure [6] shows a set of 
points for which we should find a line that supports the rest of the points, from below, and has the highest y— intersect. 
Surely, we can discard of the points on the left of the lowest point, because these points are either interior, or support 
segments of CH (DP) which have negative inclinations. Also, points on the right of the lowest point, which define lines 
with inclination greater than 1 can also be discarded. 

The last two statements are true in 2D, but not necessarily in 3D - removing points changes the convex hull, and there 
is danger that the new convex hull will introduce bogus solutions to the problem. We will see later that this is not the 
case, and removing the points is indeed safe. 



The main result is detailed in Theorem 12 which is based on the following propositions. 



Definition 7. A point p = (p x ,Py,Pz) is behind another point q = (q x , q y ,q z ) if Vx < Qxi Py < 1y an d p z > q z . 

Proposition 8. Let p mm G DP be a point with a minimal z— coordinate, and let p G DP be point which is behind p" 
Then any plane defined by n ■ (q — p) — that: 
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Figure 6: Points in red only support segments which have a negative inclination. Blue ones only support segments with 
inclinations greater than 1. 



1. Goes through p, 

2. Supports CH (DP), that is, for any point q G DP we have n ■ (q — p) > 0, and; 

3. Has a positive z— coordinate n z 

cannot be a solution to Problem^, because its dual point has negative x and/or y coordinate. 

Proof. Apply Requirement (2) to the point q = p mm ^ for which q — p equals (+a, +/3, —7) for some positive a, /3 and 7. 
The result is the relation n x a + n y /3 > n z j, which forces either of n x or n y to be positive. Now, because (n x ,n y ,n z ) is 
a normal to a plane, the corresponding plane equation must be n z ■ z (x, y) = (—n x ) x + (—n y ) y + c for some constant c. 
Because n z is positive by Requirement (3), either of the coefficients of the plane must be negative, which means the dual 
to this plane cannot be a solution to Problem [6] □ 

We can state a similar result for points which are "too steep" with respect to p min . 

Definition 9. A point p = (p Xl p y ,p z ) is too steep with respect to another point q = (q Xl q y ,q z ) if p >~ q (larger in all 
coordinates), and 

> 1 and > 1. 

Px Qx Py Qy 

Proposition 10. Let p mm € DP be a point with a minimal z— coordinate, and let p € DP be a point which is too steep 
with respect to p mm . That is, with all coordinates larger than those of p mln such that the vector p mm —p = (—a, — f3, —7), 
a, P and 7 positive, satisfies 

a < 7 and (3 < 7. 

Then any plane defined by n ■ (q — p) = that: 

1. Goes through p, 

2. Supports CH (DP), that is, for any point q G DP we have n ■ (q — p) > 0, and; 

3. Has a positive z~ coordinate n z 

cannot be a solution to Problem^ because its dual point has its x and/or y coordinate larger than 1. 
Proof. In a similar way to the proof of Proposition [8] we obtain 

(n x , n y ,n z ) ■ (-a, -f3, -7) = (-n x ) a + (-n y ) (3 - n z ~f, 
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which implies 

n z J \ n z J 

We can normalize n such that n z will be 1, without changing its sign, because n z is positive by Requirement (3): 

/ ^ a i \ P 
[-n x ) - + (-%) - > 1. 

7 7 

Now, if both a/7 and are smaller than 1, n x and n y cannot be both larger than -1. Noting that the explicit plane 
equation which the normal vector n (with n z = 1) induces is z (x, y) — (—n x ) x + (— n y ) y + c for some c, leads to the 
conclusion that at least one of the coefficients of the plane must be larger than 1, and therefore its dual point cannot be 
a solution to the linear program. □ 

We state a simple characterization of lower convex hulls : 

Lemma 11. Any plane tt defined by n ■ (q — p) — 0, which supports CH (DP) must have a positive z— coordinate of n , 
n z . 

Proof. Since tt is part of the lower convex hull, there must be a pair of points q\ and qi such that q\ is on 7r, and 

q-2 = qi + e (0, 0, 1); therefore, the condition n ■ (qi — 92) > forces n z > 0. □ 

We are now ready to show that we can safely discard of a class of constraints. 

Theorem 12. Discarding of points which are behind p mm (as defined in Proposition^, or are too steep (as defined in 
Proposition]! 0V leads to a min-max problem which is equivalent to Problemwl 



Proof. We will prove for points which are behind p mm - the proof for points which are too steep is almost identical, and 
is left for the reader. 

Let p be a point which is behind p min . First, any face which p supports and is part of CH (DP) cannot be a solution, 
therefore removing p does not change the solution to Problem [6] (apply Lemma 11 and Proposition [8]) . However, there 
still remains a possibility that removing p creates new faces in CH (DP) which results in a bogus solution to Problem [6] 
We now show it is not the case. 

We consider the incremental convex hull construction algorithm detailed in de Berg (20001. Suppose we have A = 
CH (DP\p), and would like to construct B = CH (DP). This is done by finding all the faces of A which are visible to 
p, that is, all faces which separate p from DP\p. These faces are then removed, and the hole is filled using faces that p 
supports. 

A convex hull of a set of points P is unique; one way to see this is to recall one definition of CH (P) - the intersection 
of all convex sets that contain P. This uniqueness implies that the faces that are added to the convex hull, when we 
remove p, are the faces that would have been removed if we constructed B from A using the aforementioned algorithm. 
This characterizes the faces that are added to A when we remove p - they all define planes which separate p from DP\p. 
See Figure [7] 

Formally, these planes are defined by some normal n and a point p' ', and have n ■ (q — p') > for all q G DP\p but 
n ■ (p — p 1 ) < 0. Such a plane, when translated so that its passes through p, is defined by n ■ (q — p) = 0. We note that for 
any q £ DP\p, we have 

n-(q-p) = n ■ (q - p + p' - p') 

= n-(q-p') + (-n-(p-p')) 
> 0, 

which means the translated plane goes through p, supports CH (DP) and has a positive n z component (we assumed the 
original face belonged to CH (DP)). This means it satisfies the requirements of Proposition [8] and therefore this plane 
(translated or not) cannot be a solution to Problem [6] which concludes the proof. □ 

With regard to the computational cost of these purgings, we note that they are highly-parallelizable, in addition to 
being a single-pass over the points. 
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Figure 7: Incremental convex hull construction. Left: before adding the marked point, with visible faces emphasized. 
Right: after the addition. Note that removing the marked point from the convex hull on the right image, would result in 
the convex hull that appears in the left image. 
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